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Abstract 


The proton exchange membrane fuel cell (PEMFC) has become a promising candidate for the power source of electrical vehicles because of 
its low pollution, low noise and especially fast startup and transient responses at low temperatures. A transient, three-dimensional, non-isothermal 
and single-phase mathematical model based on computation fluid dynamics has been developed to describe the transient process and the dynamic 
characteristics of a PEMFC with a serpentine fluid channel. The effects of water phase change and heat transfer, as well as electrochemical kinetics 
and multicomponent transport on the cell performance are taken into account simultaneously in this comprehensive model. The developed model 
was employed to simulate a single laboratory-scale PEMFC with an electrode area about 20cm?. The dynamic behavior of the characteristic 
parameters such as reactant concentration, pressure loss, temperature on the membrane surface of cathode side and current density during start- 
up process were computed and are discussed in detail. Furthermore, transient responses of the fuel cell characteristics during step changes and 
sinusoidal changes in the stoichiometric flow ratio of the cathode inlet stream, cathode inlet stream humidity and cell voltage are also studied and 
analyzed and interesting undershoot/overshoot behavior of some variables was found. It was also found that the startup and transient response time 
of a PEM fuel cell is of the order of a second, which is similar to the simulation results predicted by most models. The result is an important guide 


for the optimization of PEMFC designs and dynamic operation. 
© 2006 Elsevier B.V. All rights reserved. 
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1. Introduction 


A proton exchange membrane fuel cell has many useful 
characteristics, especially a short startup time at normal tem- 
perature, which makes it a promising power source for future 
transportation, such as in electrical vehicles. General speak- 
ing, PEMFCs as power sources for transportation will often 
operate dynamically, for example during startup and stop, accel- 
eration and deceleration of electrical vehicles. Thus, explicit 
knowledge of the transient behavior of a PEMFC is becom- 
ing an important research subject, and is important not only 
for the design of the controllers but also for the optimization 
of transient cell performance during changes in some operating 
cases. 
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Hitherto, there have lots of numerical investigations on PEM- 
FCs, however, most of them are based on steady and single-phase 
mathematical models and are focused on a single channel of 
a simple cell structure such as a straight channel [1-4]. By 
contrast, dynamic models and corresponding numerical simu- 
lations are rare. Some interesting and significant contributions 
have mainly been done by the researchers and groups including 
Amphlett et al. [5] and Wohr et al. [6] who mainly studied the 
dynamic responses to fuel cell temperature. Van Bussel et al. 
[7] simulated the polarization characteristics of a cell using a 
two-dimensional, dynamic model. Based on computation fluid 
dynamics model, Um and Wang [8] tried simulating the dynamic 
response of the current density to a step change in cell volt- 
age for the first time. Yerramalla et al. [9] conducted linear 
and non-linear analysis of a PEMFC system in order to know 
the whole and complex dynamic characteristics. Pathapati et al. 
[10] solved the transient response of characteristic parameters 
to a step change in current density using SIMULINK. Wang and 
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Nomenclature 


Acv 


operation area of cell (m?) 

specific surface area of the control volume (c.v.) 
(m~!) 

mass fraction of chemical species i (kg kg~!) 
constant-pressure heat capacity (Jkg~! K7!) 
effective diffusion coefficient of species i 
(m? s~!) 

Faraday constant (96487 C mol— 9) 

flux 

latent heat of water phase change (kJ kg7!) 
average current density (A m7?) 

anode exchange current density (A m7?) 


cathode exchange current density (A m~7) 


local current density (A m~?) 

permeability of porous media (m?) 
condensation rate constant (s7!) 
vapourisation rate constant ((Pa s)7!) 
molecular weight of species i (kg mol!) 
number of electrons in electrochemical reaction 
pressure of fluid (Pa) 

hydrogen partial pressure (Pa) 

oxygen partial pressure (Pa) 

stoitriometry coefficient in electrochemical reac- 
tion 

universal gas constant (8.314 J mol`! K7!) 
cell resistance 

saturation of liquid water 

source term 

entropy change for cell reaction (J K7!) 
time (s) 

temperature (K) 

intrinsical fluid velocity vector (m s7!) 
cell open-circuit voltage (V) 

mole fraction of species k 


Greek letters 
a(x, y) net water flux per proton flux 


E porosity of porous media 

n(x, y) local activation overpotential (V) 

aef effective heat conductivity 

u" effective dynamic viscosity (kg s m~7) 
0 stoichiometry flow ratio 

p density of the mixture (kg m~?) 

Omem local ionic conductivity of membrane ((Q m)~!) 
Subscripts 

a anode 

€ cathode 

i species 

k anode or cathode 

m mass 

T temperature 

u momentum 


wl liquid water 
wv water vapor 
Superscripts 

eff effective 

sat saturation 


Wang [11] estimated various time constants for important trans- 
port phenomena and analyzed characteristic parameter variation 
with time during step changes in cell voltage and the humidity of 
the inlet gas stream. Van Zee and co-workers [12,13] have con- 
ducted many experimental studies and dynamic simulations of 
a PEMFC. However, there have not been comprehensive studies 
of the characteristic parameters during startup and step changes 
in operation, in particular, transient single-phase modeling on a 
complex structure cell with a serpentine fluid channel cannot be 
found in the published literature. 

The objective of this study is two-fold. One goal is to develop 
a transient, three-dimensional, single-phase model based on 
computational fluid dynamics. The second goal with practical 
importance is to investigate the transient process of a PEMFC, 
such as to numerically analyze the transient response to fluid 
flow, reactant concentration, temperature and current density and 
so on, in a single fuel cell during startup and step changes and 
sinusoidal changes in operating cases. A complex cell structure 
with serpentine flow fields was used as the research object in 
this paper. 


2. Description of the mathematical model 


The computational domain for this model is a whole single 
cell with electrode area 4.1 cm x 5.2 cm and single serpentine 
fluid channel having 20 flow paths, as illustrated in Fig. 1. The 
bipolar plate of the fuel cell is shown in Fig. 1(a), in which 
flow paths for the reactants and products are carved. The other 
components across the membrane excluding the bipolar plate are 
shown in Fig. 1(b), including seven components: fluid channels, 
diffusion layers and catalyst layers of the anode and the cathode, 
and proton exchange membrane in the middle of the cell. 

The following simplifying assumptions have been used in 
this model development: 


(1) The flow in fuel cell is laminar everywhere. This is reason- 
able for the low velocity and low Reynolds number. 

(2) The proton membrane is impermeable to gas species. 

(3) Isothermal boundary conditions are used for the external 
wall and inlet stream. 

(4) The porous media such as the diffusion layer and the cata- 
lyst layer are isotropic and homogeneous, characterized by 
effective permeability and uniform porosity. 

(5) Transients in water accumulation in the membrane are 
neglected in this model, but are considered in Wang and 
Wang [14]. 
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(a) 


anode outlet 


cathode outlet 


(b) anode inlet 


(c) 


Fig. 1. Schematic of channel structure, components and grid partition of a sin- 
gle fuel cell. (a) Bipolar plate and channel structure, (b) components across 
membrane and (c) Numerical grid of the bipolar plate. 


(6) Water is generated in the catalyst layer in the form of vapor; 
condensation occurs after vapor partial pressure exceeds sat- 
uration vapor pressure. This is the weakest flaw in this non- 
isothermal model, and will be improved in our future work. 


2.1. Governing equation 


Numerical simulation is based on a transient, three- 
dimensional, multicomponent, non-isothermal and single-phase 


transportation model, which is acomprehensive model including 
mass, momentum, energy and chemical component transporta- 
tion. The conservation equations appear as the same form for the 
fluid channels, the diffusion layers, the catalyst layers and the 
membrane, but with different coefficients and sources because 
the catalyst layers and diffusion layers are porous media. Thus, 
the mathematical model can be described as following: 


Continuum equation 


a(ep) 
ðt 


+ V - (epu) = Sm (1) 
Momentum equation 


d(epu) 
ot 

Species conservation equation 

d(€eCi) 
ot 


+ V - (epuu) = —eVp + V- (eu! Vu) + Su (2) 


+ V- (eouC;) = V -(—pD"VCj) + S; (3) 


Energy conservation equation 


(pC pT) 


a — +V- (puCpT) = V: AFVT)+ Sr (4) 


The mass source term in continuum equation is described as 
following [15,16]: 


0, others 
Sm = XSi, catalyst layer (5) 
i 


where S; is source term of species i, such as hydrogen, water 
vapor, liquid water, oxygen. In the channels and diffusion lay- 
ers, mass source term is zero because that water vapor decreases 
equals to liquid water increases for condensation, so the total 
source is zero. Because chemical reaction can cause con- 
sumption and production of species, mass source term can be 
described as 


r 
Si = — I(x, y)MiAcy (6) 
nF 


Considering water production in the cathode catalyst layer, 
mass transfer of water through membrane for the electroosmosic 
with proton, and phase change between water vapor and liquid 
water, source term of water vapor can be written as 


1+ 2a(x, y) 
2F 


Sw,c = 


I(x, y)Mypo Acv + Swy (7) 
p 


_ ax, y) 
F 


Swa = 


I(x, y)My,oAcy + Swy, (8) 
where the source term of water vapor induced by phase change 
can be defined as 


£- (1 — s)Mp,o : CHO, sat p 
RT i 


— pX") sgn(přť — pwy) (9) 


sat 


Pwy] sgn(pwy — Pwy 


Swv, = ke 


+ kyépis(Pwy 
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Correspondingly, the source term for liquid water can be deter- 
mined as 


Swi == — sw, (10) 


p 


For the momentum equation, the source term read as follows: 


0, fluid channel 


Sy = eff (11) 
g diffusion layer/catalyst layer 


The non-zero term in above equation indicates the interaction 
force between fluid and solid base in the porous media. 

The source term for energy conservation equation can be 
written as 


Swi, x hfg, 


ST 
I(x, y? 
Omem 


For the channel and diffusion layer, latent heat of water phase 
change is considered and ohm heat of current is neglected 
because the electronic conductivity is very large; for the catalyst 
layer, the reversible and irreversible reaction heat as well as latent 
heat of water phase change are considered; for the membrane, 
ohm heating of current is considered because the resistance of 
membrane is large. 


2.2. Relation between voltage and current density 


Taking into account the voltage loss induced by ohm resis- 
tance and activation polarization, the actual cell voltage of a fuel 
cell can be explained as 


Veen = Voc — (x, y) — Ix, y) + Ric (13) 


The local activation overpotential n(x, y) is composed of 
two parts: anode activation overpotential and cathode activation 
overpotential, and the latter occupies main part, can be given by 
following equation: 


RT(x, y) I(x, y) + p(x, y) 
n(x, y) = In 
0.5F Ioo, PO, (%, y) 
RT I . 
(x, y) In (x, y): p(x, y) (14) 
1.0F Ion, PH (x, y) 


The ohm resistance includes electrode and membrane resis- 
tance, read as 


Omem Oa Oc 
Rte = +—+ (15) 
tmem tla fblc 


2.3. Initial and boundary conditions 


For species concentration and temperature, the same values 
at inlet are applied as the initial conditions. The initial condition 
of velocity can be obtained by solving flow equation without 


AS 
I(x, y) X Acy X (= + n) + Swi, X hfg, catalyst layer 
n 


reaction to reach stable state. Flow rate required for fuel and 
air electrode is a function of the desired current density and the 
stoichiometry flow ratio, which can be determined by following 
equation: 


IgA RT 


Gr=9 
k=? JF pX: 


(16) 


Furthermore, fluid velocity required for the inlet stream can be 
obtained. 

Velocity, species concentration and temperature at the inlet 
are specified according to the cell operating conditions, viz. first 


fluid channel and diffusion layer 


(12) 


: membrane 


kind boundary condition. At the outlet, pressure outlet boundary 
(assuming pressure to be zero) is employed; the other variables 
as temperature, species concentration and fluid velocity are set 
to be fully developed, that is, the gradient of these variables is 
equal to zero. For the interface between fluid domain and solid 
domain, the flux of species is equal to zero, no-slip applies to 
velocity condition and the coupled boundary condition is set 
for temperature. The temperature on the external surface of the 
bipolar plate is specified according to cell operation environment 
temperature, such as 353.15 K used in this paper. 


2.4. Numerical procedure 


The whole solution includes three steps corresponding to 
three main processes of cell operation. First of all, the com- 
putation without reaction is conducted until the solution is 
convergent, which corresponds to the pre-feed reactants and pre- 
heating process for fuel cell operation. Secondly, the transient 
startup process computation using results of previous process as 
initial value is conducted, and until the new steady-state value is 
obtained. The reaction source terms are added for this compu- 
tation, which corresponds to the cell startup process after being 
connected to work. Finally, computations of the dynamic process 
during step and sinusoidalchanges in operating conditions are 
conducted. The purpose of this computation is to know the tran- 
sient response characteristics of fuel cell, and analyze dynamic 
performance of fuel cell. 

All the coupled, non-linear model equations, including mass, 
momentum, species, temperature conservation equations, were 
converted to the finite-difference form by the control volume 
method. Stringent numerical tests are performed to ensure that 
the solution is independent of the grid size. As a result, the 
82 x 20 x 104 mesh can provide sufficient spatial resolution. 
The uniform grid is partitioned in direction y and z, but non- 
uniform in direction x. Fig. 1(c) shows the computation grid that 
was used for the bipolar plate. The solutions of the velocity com- 
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ponent equations are obtained in a staggered control volume. 
The coupled fluid velocity fields and pressure are solved by 
semi-implicit method for pressure-linked equations consistent 
(SIMPLEC) arithmetic. Then, the species and energy trans- 
port equations are computed after fluid velocities are obtained. 
Updating parameters based the computed results; the external 
iteration is conducted until the resolution is convergent. In this 
paper, the convergent criterion is that relative errors of all vari- 
ables are less than 1074. In order to obtain the polarization 
curves, the method employed in this paper is to solve average cur- 
rent for a given voltage, which can guarantee solution stability 
and convergence. The time step for transient solution is 0.02 s. 

This model is implemented into the commercial CFD code 
Fluent 6.1, with custom developed user-subroutines that con- 
sider the physicochemical processes associated with PEM fuel 
cells. 


2.5. Physical properties 


Mass diffusion coefficient of species i can be obtained 
through mass weight average of binary diffusion coefficient: 


Di = Simi : Di, (17) 
a 


where the binary diffusion coefficient D; j is defined as [17]: 


T 
Di; = 3.64 x 1078 | ——— |- (Toi - Te p? 
PËi,j AA (Tei + Te, j) 
i 1 1\ 12 
‘(De,i © Pe, j) B, (z + iz) (18) 
i j 


And for the porous media as diffusion layer, the diffusion 
coefficient should be correlated by considering the effect of 
porosity and tortuosity factor. Furthermore, liquid saturation 
is considered to affect species diffusion coefficient. Thus, 
the effective mass diffusion coefficient can be determined 
as: 


eff _ (1 —s)- Di, channel 


a (1—s)D;-e!°, diffusion layer and catalyst layer 


(19) 


The effective heat conductivity can be defined as [18] 


kg, channel 
=i 
Aett = 4 —2ks + ase + a *) , porous media 
64h She 
ks, solid plate 


(20) 


Local membrane conductivity Omem not only affects current 
density of the cell, but also determines temperature distribution 
inside the cell. Thus, it is important to obtain an accurate pre- 
diction of membrane conductivity, which is a function of water 
content and temperature of membrane. We adapt a widely used 


experiential equation summarized by Springer and Zawodinski 


[2]: 


ref 0.005139A — 0.00326, A> 1 
o = (21) 
constant, à<l1l 


where am is the reference protonic conductivity under temper- 
ature 303 K and A is the membrane water content. It can be 
correlated by considering temperature effect: 


“ap [1288 (335-7) 
oT) = ot exp | 1268 | — — — (22) 


303 T 
The relation between water content A at membrane surface and 
water vapor activity can be described by experiential equation 
as [19]: 


0.043 + 17.81a — 39.85a? + 36.00, O<a<1 
144+ 1.4(a — 1), l<a<3 
A= (23) 
16.8, a>3 
22, liquid water existence 


where a stands for water vapor activity at membrane surface. 
Water content is assumed to increase linearly from 14 under 
saturation to 16.8 under activity 3, because the experimental 
data for membrane content under activities 1-3 is scanty. Water 
vapor activity at membrane surface a can be defined as: 

ea (24) 


ps 


where p* is saturation vapor pressure corresponding to local 
temperature. The experiential expression equation between sat- 
uration vapor pressure and temperature is obtained by fitting 


Table 1 

Physics parameters and base conditions 

Quantity Value 
Operating temperature (K) 353.15 
Length of the channel (cm) 5.0 
Channel width (cm) 0.1 
Width of plate out channel (cm) 0.1 
Channel thickness (cm) 0.1 
Bipolar plate width (cm) 41 x 0.1 
Diffusion layer thickness (cm) 0.025 
Catalyst layer thickness (cm) 0.0028 
Membrane thickness (cm) 0.005 
Shoulder width (cm) 0.1 
Operating voltage (V) 0.65 
Pressure of anode inlet (kPa) 1 x 101.3 
Pressure of cathode inlet (kPa) 1 x 101.3 
Exchange current density of anode (A m7?) 100 
Exchange current density of cathode (A m~?) 1000 
Specific area of catalyst layer (m7!) 1.4x 10° 
Velocity at the anode inlet (m s7!) 2.46 
Velocity at the cathode inlet (ms~!) 13.15 
Humidity of anode inlet stream (%) 80 
Humidity of cathode inlet stream (%) 100 
Condensation rate (s~!) 1.0 
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saturation vapor pressure under different temperature [3]: 


logiop* = —2.1794 + 0.02953T — 9.1837 x 10757? 


+ 1.4454 x 10-7 7F? (25) 


3. Results and discussion 


Parameters of the fuel cell dimensions and the base operating 
case used in the computation of this study are listed in Table 1. 
Other physical parameters can be obtained from Bernardi 
and Verbrugge [1], or solved by method introduced in this 
literature. 
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3.1. Dynamic characteristics during startup process 


The variations of hydrogen and oxygen concentration in the 
x-z plane at half height of channel with time during startup 
process are illustrated in Figs. 2 and 3. First of all, hydrogen con- 
centration increases and oxygen concentration decreases along 
with flow for the electrochemical reaction. Although, both of 
them are consumed in the reaction, but there exist net water trans- 
porting from anode to cathode along with proton transfer, which 
cause the hydrogen concentration increase relatively because 
the molecular weight of water is much larger than hydrogen. 
Comparing these two pictures, we can also find that the time 
for hydrogen approaching steady-state value is longer, about 


0.05 : 
0.04 
0.240 
0.234 
0.228 
0.222 
0.03 0.216 
g 0.210 
N 0.204 
N 0.199 
0.193 
0.02 0.187 
0.181 
0.175 
0.169 
0.01 0.163 
0.157 
0 pu T L s UN u 
0 0.01 0.02 0.03 0.04 
(b) x/m 
0.05 
0.04 0.258 
0.251 
0.244 
0.236 
0.03 0.229 
f=] 0.222 
N] 0.215 
N 0.207 
0.200 
0.02 0.193 
0.186 
0.179 
0.171 
0.01 0.164 
0.157 
0 
0 0.01 0.02 0.03 0.04 
(d) x/m 
0.05 
0:04 0.260 
0.253 
0.245 
.23 
0.03 0231 
g 0.223 
0.216 
N 0.208 
0.201 
0.02 0.194 
0.186 
0.179 
0.172 
0.01 0.164 
0.157 
0 
0 0.01 0.02 0.03 0.04 
(f) x/m 


Fig. 2. Variations of hydrogen concentration in the x-z plane at half height of anode channel with time (a) t=0.3 s, (b) f=0.6s, (c) t=0.9s, (d) t= 1.25, (e) f=1.5s 


and (f) t=1.8s. 
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Fig. 3. Variations of oxygen concentration in the x-z plane at half height of cathode channel with time (a) t=0.3 s, (b) t=0.6s, (c) t=0.9s, (d) t=1.2s, (e) t=1.5s 


and (f) t=1.8s. 


1.2-1.5 s, in contrast, the time for oxygen approaching steady- 
state value is relatively shorter, about 0.6—0.9 s. The reason for 
this is mainly that the velocity in the cathode channel is faster 
than that in the anode channel, so the reactants can be provided 
by convection more quickly, i.e., the convection mass transfer 
takes up a major position, although the diffusion coefficient of 
hydrogen is larger than that of oxygen. We can also say that 
mass transfer is more important than the reaction rate because 
the oxygen reaction rate is much slower than the hydrogen reac- 
tion, i.e., the time constant of electrochemical reaction is very 
short in contrast to that of mass transfer. 

Fig. 4 plots variations of the local current density with 
time during the startup process. Similarly, local current density 


decreases from the inlet to the outlet during the consumption 
of reactants. The time in approaching the steady-state value is 
also very short, about 0.6-0.9 s. This time is very close to the 
time required by oxygen concentration approaching steady state, 
which indicates that the time constant of the cell is determined by 
oxygen transfer. Because the current density is the main reflec- 
tion of cell operation characteristics, it enough to imply that the 
startup time of the fuel cell is very short and of the order of 
second. Thus, a PEMFC can match the requirement as a power 
sources for transportation such as an electric vehicle. Actually, 
for a whole powering system, there exist many secondary sys- 
tems each requiring some response time and the startup time of a 
PEMFC system is a second order effect; however, this is beyond 
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Fig. 4. Variations of local current density with time (a) t=0.3 s and (b) t=0.6s, (c) t=0.9s, (d) t= 1.2, (e) t= 1.5 and (£) t=1.8s. 


the research range of this paper, but it is an important subject for 
our future research. 

Fast startup under normal temperature is the prominent 
advantage of a PEMFC, which is different from other fuel cells. 
The computed variations of cell characteristics with time dur- 
ing the startup process are shown in Fig. 5: (a) pressure loss of 
the anode, (b) pressure loss of the cathode, (c) average temper- 
ature on the membrane surface of the cathode side, (d) oxygen 
concentration in the catalyst layer of the cathode, (e) hydrogen 
concentration in the catalyst layer of the anode and (f) average 
current density. Pressure loss for the anode is smaller than for 
the cathode; the main reason is that the inlet velocity of the cath- 
ode gas stream is higher than that of the anode. Furthermore, the 
pressure losses of both anode and cathode are smaller than with- 


out reaction. The average temperature on the membrane surface 
of the cathode side increases with time until a stable state is 
reached; this is because the reaction emits thermal energy. Sim- 
ilar to Fig. 2, the hydrogen mass fraction in the anode increases 
with time, this is because more water moves from the anode 
to the cathode during electrochemical reaction, then causes the 
mass fraction of hydrogen to increase relatively. The oxygen 
mass fraction decreases during consumption, which is similar to 
the decrease of local current density. Also, it can be found that 
there exists an interesting undershoot phenomenon for pressure 
loss of the anode and oxygen concentration of the cathode, but 
other variables have not shown this phenomenon. Altogether, 
the characteristic parameters reach a steady-state value in less 
than 2s, i.e., the startup time of PEMFC is very short. 
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Fig. 5. Characteristic parameters variation during startup process. (a) Anode pressure loss, (b) cathode pressure loss, (c) average temperature on the membrane 
surface of the cathode side, (d) oxygen concentration in the catalyst layer of cathode side, (e) hydrogen concentration in the catlyst layer of anode side and (f) average 


current density (A m7). 
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3.2. Dynamic response during step change in operating 
cases 


The following results were computed with a step change in 
one operating parameter, keeping the others constant. The tem- 
perature on the membrane surface of the cathode side was higher 
than in the other domains, so it was used as the research object 
in the following discussion. 

Fig. 6 illustrates the dynamic response of the characteristic 
parameters during a step change in the stoichiometry flow ratio 
of the cathode inlet, which changes between 1.0 and 1.5 period- 
ically. Fig. 6(a) shows a step change in the stoichiometry flow 
ratio, and Fig. 6(b and c) are the dynamic responses of the cor- 
responding temperature of the membrane surface of the cathode 
side and the average current density of the fuel cell. It was found 
that the change in stoichiometry flow ratio had little effect on 
the temperature, about 0.02 K. The main reason is that the sto- 
ichiometry flow ratio increases, current density increases, and 
more heat was be produced; however, on the other hand more 
thermal energy can be removed by the higher velocity gas stream, 
serving as an effective refrigerant, and thus the total increase 
of temperature is very small. Increasing the inlet stoichiome- 
try flow ratio can increase the cell performance effectively, i.e., 
the change of current density is larger, 5688-5768 A m~7. The 
reason is that the reactants can be provided in time and the 
products can be removed faster. Response times of the cell to 
a step change between both stoichiometry flow ratios are very 
close, about 0.5 s. This indicates that the fuel cell can stabilize 
in a short time during a step change in the stoichiometric flow 
ratio. There exists a small overshoot of the average current den- 
sity when the stoichiometry flow ratio changes from large to 
small, by contrast, an undershoot of the average current density 
occurs when the stoichiometry flow ratio changes from small to 
large. The temperature increases when the stoichiometric flow 
ratio increases, and vice versa, however, there does not exist 
any overshoot and undershoot phenomenon during the change 
process. 

Fig. 7 shows the dynamic response behavior during a step 
change in the cathode inlet stream humidity, i.e., the rela- 
tive humidity changes between 100 and 80% in the cycle. 
Fig. 7(a) plots the step change in the relative humidity, Fig. 7(b 
and c) show the corresponding dynamic responses of the 
temperature on the membrane surface of the cathode side 
and the average current density. It was found that the tem- 
perature changes about 0.1 K during the step change in the 
relative humidity, however the current density had relatively 
large change, from 5688 to 6400 Am~?. When the humid- 
ity changes from large to small, the average current density 
increased and there existed a weak overshoot phenomenon. 
When the humidity changed from small to large, the aver- 
age current density decreased and there existed an undershoot 
behavior. The temperature increased when humidity decreased, 
and vice versa, however, there did not exist any overshoot or 
undershoot phenomenon. The temperature will achieve a steady 
state in a shorter time than the current density during the step 
change in the cathode inlet stream humidity, of 0.7 and 1.4s, 
respectively. 
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Fig. 6. Dynamic response of fuel cell during step change in stoichiometry flow 
ratio at cathode inlet. (a) Stiochiometry flow ratio at cathode inlet, (b) average 
temperature on the membrane surface of the cathode side and (c) average current 
density (A m~). 


Fig. 8 shows the dynamic response during a step change in 
the operating voltage, which changes between 0.65 and 0.7 V in 
the cycle. Fig. 8(a) plots the step change in the operating volt- 
age, and Fig. 8(b and c) illustrate the corresponding dynamic 
response of the temperature on the membrane surface of the 
cathode side and the average current density. It can be seen from 
the figures that the temperature has a small change of about 
0.1 K during the step change in the operating voltage, however 
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Fig. 7. Dynamic response of fuel cell during step change in cathode inlet stream 
humidity. (a) Cathode inlet stream humidity, (b) average temperature on the 
membrane surface of the cathode side and (c) average current density (A m7”). 


the current density change was relatively large, from 4300 to 
5688 A m™?. In the higher cell voltage case, the current density 
was lower and the temperature was lower for small ohmic heat- 
ing effect. During a step change in cell voltage, there existed an 
obvious overshoot and undershoot phenomenon of current den- 
sity, of about 200 A m2; however, the temperature did not show 
this phenomenon. The current density achieved a steady state in 
a shorter time, about 0.5 s, however the temperature achieved 
steady state in about 1.8 s. 
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Fig. 8. Dynamic response of fuel cell during step change in operating voltage. 
(a) Operating voltage, (b) average temperature on the membrane surface of the 
cathode side and (c) average current density (A m7). 


3.3. Dynamic response behavior during sinusoidal change 
in operating cases 


The following results are computed with a sinusoidal change 
in one characteristic parameter, keeping others constant. The 
purpose of this computation is to know the effect of a sinusoidal 
change on the cell performance. 
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Fig. 9. Dynamic response of fuel cell during sinusoidal change in stoichiometry 
flow ratio at cathode inlet. (a) Stoichiometry flow ratio of cathode inlet, (b) 
average temperature on the membrane surface of the cathode side and (c) average 


current density. 
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Fig. 10. Dynamic response of fuel cell during sinusoidal change in cathode inlet 
stream humidity. (a) cathode inlet stream humidity, (b) average temperature on 
the membrane surface of the cathode side and (c) average current density. 
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Fig. 11. Dynamic response of fuel cell during sinusoidal change in operating 
voltage. (a) Operating voltage, (b) average temperature on the membrane surface 
of the cathode side and (c) average current density. 


Fig. 9 illustrates the dynamic response during a sinu- 
soidalchange in the stoichiometry flow ratio of cathode inlet 
stream. Fig. 9(a) shows a sinusoidal change of the stoichiom- 
etry flow ratio of the cathode inlet stream, and Fig. 9(b and c) 
are the dynamic response behaviors of the corresponding tem- 
perature on the membrane surface of the cathode side and the 
average current density. It can be seen that the temperature and 
current density increase along with increase of the stoichiome- 
try flow ratio, vice versa, decreases along with the decrease of 
the stoichiometry flow ratio. The process from high stoichiom- 
etry flow ratio to low stoichiometry flow ratio and from low 
stoichiometry flow ratio to high stoichiometry flow ratio is not 
very symmetrical; the change of value close to the upper limit 
is slower than that close to lower limit. Comparing to a step 
change, the upper limit and lower limit of temperature and cur- 
rent density are very close, but there does not exist any overshoot 
and undershoot phenomenon for the current density for the slow 
change of stoichiometry flow ratio. 

Fig. 10 shows the dynamic response behavior during a sinu- 
soidal change in the cathode inlet stream humidity, i.e., relative 
humidity changes from 100 to 80% in cycle. Fig. 10(a) plots 
the sinusoidal change of the relative humidity, Fig. 10(b and c) 
are the corresponding dynamic responses of the temperature on 
the membrane surface of the cathode side and the average cur- 
rent density. It was found that the temperature changed little, 
about 0.1 K, however the current density had a relatively large 
change, from 5688 to 6400 A m~*. When the humidity changed 
from high to low, the average current density increased, and the 
temperature increased, and vice versa, the current density and 
temperature decreased. Comparing to a sinusoidal change in the 
stoichiometry flow ratio, the process from high humidity to low 
humidity and from low humidity to high humidity was almost 
symmetrical. 

Fig. 11 shows dynamic response during sinusoidal change 
in cell voltage, between 0.65 and 0.7 V in the cycle. Fig. 11(a) 
plots the sinusoidal change in cell voltage, and Fig. 11(b and 
c) illustrate the corresponding dynamic response of the tem- 
perature on the membrane surface of the cathode side and the 
average current density. It was found that the temperature was lit- 
tle changed, less than 0.1 K, however the current density changed 
relatively much more, from 4300 to 5688 A m~?. Under a higher 
cell voltage, the current density is lower and the temperature was 
lower for a small ohmic heating effect. This is a suitable polar- 
ization relationship between the current density and operation 
voltage. 


4. Conclusions 


A transient, three-dimensional and single-phase mathemat- 
ical model was developed to investigate the dynamic response 
behavior of the characteristic PEMFC parameters such as tem- 
perature and current density during the startup process and step 
change and sinusoidal changes in stoichiometric flow ratio of the 
cathode inlet stream, the cathode inlet stream humidity and the 
cell voltage. Interesting undershoot/overshoot behavior of some 
of the variables was found. The computed results indicate that 
the dynamic response processes of PEMFCs are very fast, only 
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about 1-2 s, which is similar to most other reported simulation 
results, however this is significantly different from the time con- 
stant in Wang and Wang [11], the main reason is that the time 
change of liquid water in the membrane is considered in their 
paper. The results also show that this kind of low temperature 
fuel cell can be suitable for equipment requiring a quick startup 
power source such as an electrical vehicle. 


Acknowledgements 


The authors wish to express their thanks to Natural Science 
Foundation of Zhejiang Province of P.R. China (Grant No. 
Y106161). 


References 


[1] D.M. Bernardi, M.W. Verbrugge, A mathematical model of the solid- 
polymer-electrolyte fuel cell, J. Electrochem. Soc. 139 (1992) 2477- 
2491. 

[2] T.E. Springer, T.A. Zawodinski, Polymer electrolyte fuel cell model, J. 

Electrochem. Soc. 138 (1991) 2334-2342. 

[3] S. Dutta, S. Shimpalee, J.W. Van Zee, Three-dimensional numerical simu- 

lation of straight channel PEM fuel cells, J. Appl. Electrochem. 30 (2000) 

135-146. 

[4] G.L. Hu, J.R. Fan, S. Chen, Y.J. Liu, K.F. Cen, Three-dimensional numer- 
ical analysis of proton exchange membrane fuel cells (PEMFCs) with 
conventional and interdigitated flow fields, J. Power Sources 136 (1) (2004) 
1-9. 

[5] J.C. Amphlett, et al., A model predicting transient responses of proton 
exchange membrane fuel cells, J. Power Sources 61 (1996) 183. 


[6] M. Wohr, Bolwin, W. Schnurnberger, Dynamic modeling and simulation 
of a polymer membrane fuel cell including mass transport limitations, Int. 
J. Hydrogen Energy 23 (3) (1998) 213-218. 

[7] H.P.L.H. Van Bussel, et al., Dynamic model of solid polymer fuel cell water 

management, J. Power Sources 71 (1998) 218-222. 

[8] S. Um, C.Y. Wang, Computational fluid dynamics modeling of proton 

exchange membrane fuel cells, J. Electrochem. Soc. 147 (12) (2000) 

4485-4493. 

[9] S. Yerramalla, A. Davari, A. Faliachi, et al., Modeling and simulation of the 
dynamic behavior of a polymer electrolyte membrane fuel cell, J. Power 
Sources 124 (2003) 104-113. 

[10] P.R. Pathapati, et al., A new dynamic model for predicting transient phe- 

nomena in a PEM fuel cell system, Renewable Energy 30 (2005) 1-22. 

[11] Y. Wang, C.Y. Wang, Transient analysis of polymer electrolyte fuel cell, 
Electrochim. Acta 50 (2005) 1307-1315. 

[12] S. Kim, S. Shimpalee, J.W. Van Zee, The effect of reservoirs and fuel 
dilution on the dynamic behavior of a PEMFC, J. Power Sources 137 (2004) 
43-52. 

[13] S. Kim, S. Shimpalee, J.W. Van Zee, The effect of stoichiometry on dynamic 
behavior of a proton exchange membrane fuel cell (PEMFC) during load 
change, J. Power Sources 135 (2004) 110-121. 

[14] Y. Wang, C.Y. Wang, Dynamics of polymer electrolyte fuel cells undergo- 
ing load changes, Electrochim. Acta 51 (2006) 3924-3933. 

[15] C.Y. Wang, Fundamental models for fuel cell engineering, Chem. Rev. 104 
(2004) 4727-4766. 

[16] Y. Wang, C.Y. Wang, Modeling polymer electrolyte fuel cells with large 
density and velocity changes, J. Electrochem. Soc. 152 (2005) A445-A453. 

[17] R. Bird, et al., Transport Phenomena, Wiley, New York, 1960. 

[18] V. Gurau, H. Liu, S. Kakac, Two-dimensional model for proton exchange 
membrane fuel cells, AICHE J. 44 (1998) 2410-2422. 

[19] T.A. Zawodzioski, C. Derouin, S. Radzinski, et al., Water uptake by and 

transport through Nafion 117 mambrane, J. Electrochem. Soc. 140 (1993) 

1041. 


